#include<iostream>
#include<stdlib.h>
#include<iomanip>
#include<ctime>
#include<windows.h>
#include <xmmintrin.h> //SSE
#include <emmintrin.h> //SSE2
#include <pmmintrin.h> //SSE3
#include <tmmintrin.h> //SSSE3
#include <smmintrin.h> //SSE4.1
#include <nmmintrin.h> //SSSE4.2
#include <immintrin.h> //AVX、AVX2
using namespace std;
const int n = 500;
void init(float A[n][n])//n*n矩阵赋值
{
    srand((unsigned)time(NULL));
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < i; j++)
        {
            A[i][j] = 0;
        }
        A[i][i] = 1.0;
        for (int j = i + 1; j < n; j++)
        {
            A[i][j] = rand() % 10000;//模取一万保证矩阵数据不会过大
        }
    }
    for (int k = 0; k < n; k++)
    {
        for (int i = k + 1; i < n; i++)
        {
            for (int j = 0; j < n; j++)
            {
                A[i][j] += A[k][j];
                A[i][j] = float(int(A[i][j]) % 10000);//同理
            }
        }
    }
}

void display(float A[n][n])
{
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            cout << setw(11) << A[i][j] << " ";
        }
        cout << endl;
    }
}

void Gauss_normal(float** A)//普通高斯消元算法
{
    for (int k = 0; k < n; k++)
    {
        for (int j = k + 1; j < n; j++)
        {
            A[k][j] = A[k][j] / A[k][k];
        }
        A[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            for (int j = k + 1; j < n; j++)
            {
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            }
            A[i][k] = 0;
        }
    }
}
void clone(float A[n][n], float** B)
{
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            B[i][j] = A[i][j];
        }
    }
}
void Gauss_normal_sse(float** A)//高斯消元sse优化算法
{
    for (int k = 0; k < n; k++)
    {
        __m128 vt = _mm_set1_ps(A[k][k]);
        int j;
        for (j = k + 1; j + 4 <= n; j += 4)
        {
            __m128 va = _mm_loadu_ps(&A[k][j]);
            va = _mm_div_ps(va, vt);
            _mm_storeu_ps(&A[k][j], va);
        }
        for (; j < n; j++)
        {
            A[k][j] = A[k][j] / A[k][k];
        }
        A[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            __m128 vaik = _mm_set1_ps(A[i][k]);
            for (j = k + 1; j + 4 <= n; j += 4)
            {
                __m128 vakj = _mm_loadu_ps(&A[k][j]);
                __m128 vaij = _mm_loadu_ps(&A[i][j]);
                __m128 vx = _mm_mul_ps(vakj, vaik);
                vaij = _mm_sub_ps(vaij, vx);
                _mm_storeu_ps(&A[i][j], vaij);
            }
            for (; j < n; j++)
            {
                A[i][j] = A[i][j] - A[k][j] * A[i][k];
            }
            A[i][k] = 0;
        }
    }
}
void Gauss_normal_avx(float** A)//高斯消元avx优化算法
{
    for (int k = 0; k < n; k++)
    {
        __m256 vt = _mm256_set1_ps(A[k][k]);
        int j;
        for (j = k + 1; j + 8 <= n; j += 8)
        {
            __m256 va = _mm256_loadu_ps(&A[k][j]);
            va = _mm256_div_ps(va, vt);
            _mm256_store_ps(&A[k][j], va);
        }
        for (; j < n; j++)
        {
            A[k][j] = A[k][j] / A[k][k];
        }
        A[k][k] = 1.0;
        for (int i = k + 1; i < n; i++)
        {
            __m256 vaik = _mm256_set1_ps(A[i][k]);
            for (j = k + 1; j + 8 <= n; j += 8)
            {
                __m256 vakj = _mm256_loadu_ps(&A[k][j]);
                __m256 vaij = _mm256_loadu_ps(&A[i][j]);
                __m256 vx = _mm256_mul_ps(vakj, vaik);
                vaij = _mm256_sub_ps(vaij, vx);
                _mm256_store_ps(&A[i][j], vaij);
            }
            for (; j < n; j++)
            {
                A[i][j] = A[i][j] - A[k][j] * A[i][k];
            }
            A[i][k] = 0;
        }
    }
}
int main()
{
    float A[n][n];
    init(A);
    double time1_sum = 0;
    int count1 = 0;
    for (int l = 0; l < 100; l++)//测量平凡算法时间，运行十次求平均值
    {
        double time1 = 0;
        float** B = new float* [n];
        for (int i = 0; i < n; i++)
        {
            B[i] = new float[n];
        }
        clone(A, B);
        long long freq, head, tail;
        QueryPerformanceFrequency((LARGE_INTEGER*)&freq);
        QueryPerformanceCounter((LARGE_INTEGER*)&head);
        Gauss_normal(B);
        QueryPerformanceCounter((LARGE_INTEGER*)&tail);
        time1 = (tail - head) * 1000 / freq;
        if (time1 != 0)count1++;
        time1_sum += time1;
    }
    cout << "普通算法总运行时间：" << time1_sum / count1 << "ms" << endl;

    double time2_sum = 0;
    int count2 = 0;
    for (int l = 0; l < 100; l++)//测量sse优化算法时间
    {
        double time2 = 0;
        float** B = new float* [n];
        for (int i = 0; i < n; i++)
        {
            B[i] = new float[n];
        }
        clone(A, B);
        long long freq, head, tail;
        QueryPerformanceFrequency((LARGE_INTEGER*)&freq);
        QueryPerformanceCounter((LARGE_INTEGER*)&head);
        Gauss_normal_sse(B);
        QueryPerformanceCounter((LARGE_INTEGER*)&tail);
        time2 = (tail - head) * 1000 / freq;
        if (time2 != 0)count2++;
        time2_sum += time2;
    }
    cout << "普通算法sse总运行时间：" << time2_sum / count2 << "ms" << endl;

    double time3_sum = 0;
    int count3 = 0;
    for (int l = 0; l < 100; l++)//测量avx优化算法时间
    {
        double time3 = 0;
        float** B = new float* [n];
        for (int i = 0; i < n; i++)
        {
            B[i] = new float[n];
        }
        clone(A, B);
        long long freq, head, tail;
        QueryPerformanceFrequency((LARGE_INTEGER*)&freq);
        QueryPerformanceCounter((LARGE_INTEGER*)&head);
        Gauss_normal_sse(B);
        QueryPerformanceCounter((LARGE_INTEGER*)&tail);
        time3 = (tail - head) * 1000 / freq;
        if (time3 != 0)count3++;
        time3_sum += time3;
    }
    cout << "普通算法avx总运行时间：" << time3_sum / count3 << "ms" << endl;
}
